Understanding and leveraging phenotypic plasticity during metastasis formation

Cancer metastasis is the process of detrimental systemic spread and the primary cause of cancer-related fatalities. Successful metastasis formation requires tumor cells to be proliferative and invasive; however, cells cannot be effective at both tasks simultaneously. Tumor cells compensate for this trade-off by changing their phenotype during metastasis formation through phenotypic plasticity. Given the changing selection pressures and competitive interactions that tumor cells face, it is poorly understood how plasticity shapes the process of metastasis formation. Here, we develop an ecology-inspired mathematical model with phenotypic plasticity and resource competition between phenotypes to address this knowledge gap. We find that phenotypically plastic tumor cell populations attain a stable phenotype equilibrium that maintains tumor cell heterogeneity. Considering treatment types inspired by chemo- and immunotherapy, we highlight that plasticity can protect tumors against interventions. Turning this strength into a weakness, we corroborate current clinical practices to use plasticity as a target for adjuvant therapy. We present a parsimonious view of tumor plasticity-driven metastasis that is quantitative and experimentally testable, and thus potentially improving the mechanistic understanding of metastasis at the cell population level, and its treatment consequences.

transitions. Phenotype i grows at rate r i and transitions to the adjacent more mesenchymal-like type at rate T EM and to the adjacent more epithelial-like type at rate T M E . Resource competition is modeled with the term r i X K where X = N j=1 x j is the total population abundance, and K is the carrying capacity. The epithelial phenotype can only transition to the more mesenchymallike adjacent hybrid phenotype, and the mesenchymal phenotype can only transition to the epithelial-like adjacent hybrid phenotype; thus, the epithelial and mesenchymal phenotypes are the terminal phenotypes. Here, we show only one hybrid phenotype, but we also investigate the eect of a larger number of hybrid phenotypes.
Phenotype transitions generate and maintain heterogeneity. 134 Assuming that cells can transition into adjacent epithelial or mesenchymal phenotypes generates 135 a phenotype distribution of abundances along the epithelial-mesenchymal trait axis. Without 136 treatment, m D = m I = 0, the model (see Section Model) has two equilibria. One equilibrium 137 is the state where the tumor vanishes, x i = 0 for i = 1, 2, . . . , N . This equilibrium is unstable, 138 implying that a tumor can always progress and increase in abundance without treatment. The 139 second equilibrium is the state described by 140 x for i = 1, 2, . . . , N.
(1) 141 This equilibrium is stable and describes the coexistence of all phenotypes (see supplementary text 142 for a proof of global stability). This stability implies that phenotypic heterogeneity is generated 143 and maintained in tumor populations featuring phenotypic transitions. Patient-derived breast 144 cancer cell populations exhibit similar behavior in-vitro 34,35 . 145 The sum of all phenotype abundances at the coexistence equilibrium is K. The stable 146 distribution of phenotypes at the coexistence equilibrium depends on the transition bias λ and 147 the number of phenotypes N . Notably, the stable phenotype distribution does not depend on 148 the growth rate. Thus, only the transitions and not the growth dynamics decide the phenotype 149 abundances at equilibrium. Additionally, this growth independence shows that the particular 150 choice of the decline of growth rates from epithelial to mesenchymal phenotypes does not aect 151 the equilibrium phenotype distribution. 152 The transition bias controls the stable phenotype distribution. 153 When cells change to a more epithelial or mesenchymal phenotype with the same probability, 163 i.e., when there is no transition bias (λ = 0), all phenotypes are at equal abundance and the 164 stable distribution is a uniform distribution (Fig. 2). The mesenchymal phenotype (M) is the 165 most abundant when there is a transition bias towards the mesenchymal phenotypes (λ > 0). 166 On the other hand, for a transition bias towards the epithelial phenotypes (λ < 0), the epithelial 167 phenotype (E) is the most abundant. The phenotypic heterogeneity of the tumor population 168 is highest when there is no transition bias and λ = 0 (Fig. S1). A high variance corresponds 169 to presence of a lot of dierent phenotypes reducing the ecacy of standard cancer treatments. 170 The heterogeneity decreases and vanishes at the bias extremes λ → ±1 where transitions become 171 unidirectional. 172 The presence of phenotypes far o the average population phenotype with very low 173 abundances is captured by the third central moment of the phenotype distribution. These 174 phenotypes are essential for residual disease in some cancers, and thus an important potential 175 treatment target. We nd that the third central moment is an odd function of the transition 176 bias, vanishing at λ = 0 and at high transition bias (Fig. S1). The number of phenotypes N 177 quantitatively aects the moments of the stable phenotype distribution but does not aect their 178 shapes qualitatively (Fig. S1). We will thus x N = 3 phenotypes in the remainder of the study 179 for illustration. (1)) relative to the carrying capacity K for a xed value of the transition bias λ and the number of phenotypes N . The stable phenotype distribution changes with the transition bias λ but remains qualitatively unaected by changing the number of phenotypes N . The stable distribution is uniform when there is no transition bias to either epithelial or mesenchymal-like phenotypes, i.e., λ = 0. λ < 0 depicts a transition bias towards epithelial-like phenotypes and leads to a relative increase in epithelial cells. Conversely, λ > 0 results in a transition bias towards mesenchymal-like phenotypes and causes a relative increase in mesenchymal cells. to the primary site, the stable equilibrium is approached by a relative increase of epithelial-like 213 cells at the secondary site (Fig. 4, second column) during metastasis formation. 214 After a perturbation from the coexistence equilibrium, e.g., by treatment, the stable 215 phenotype distribution is restored, compensating the eect of the perturbation (Fig. 4, 216 columns three to ve). Hence, the phenotypically plastic tumor cell population can maintain 217 heterogeneity. Our model thus exhibits site-and treatment-specic phenotype shifts that result 218 only from the selection during dissemination or treatment. These shifts do not require specic 219 microenvironmental factors, although the microenvironment is a crucial component contributing 220 to the initial conditions that determine the direction of the shift. 221 Treatment type aects the phenotype distribution of the tumor. 222 Next, we tested how the phenotype distribution changes during and after growth-dependent 233 and growth-independent treatment (Fig. 5). Growth-dependent treatment targets proliferating 234 cells, thus epithelial cells experience higher treatment-induced mortality (Fig. 5a), and the mean 235 of the distribution shifts towards the mesenchymal phenotype (Fig. 5c). Growth-independent 236 treatment targets all phenotypes equally, not changing the phenotype distribution directly. The 237 room for growth is lled by epithelial-like cells as they grow fast (Fig. 5b)   : Plasticity-driven approach to stable phenotype distribution from dierent initial conditions. Depending on the initial condition, the approach to the stable phenotype distribution proceeds along dierent paths (top four rows, the vertical axis represents time). The mean population phenotype over time is shown in the last row with its asymptote (dashed line).  Figure 5: Growth-dependent and growth-independent treatment can transiently alter the phenotype distribution of a tumor. During treatment, the tumor burden reduces (panels a and b), and the mean phenotype of the tumor changes if the transition speed c is small but restores to the untreated mean phenotype after treatment (panels c and d). Treatment is applied between t = 0 and t = 10. Afterward, regrowth is tracked until t = 1000. The violet and green horizontal bars indicate growth-dependent and growth-independent treatments. Growthdependent treatment shifts the mean of the phenotype distribution towards the mesenchymal phenotype as it exerts higher mortality on epithelial cells. Conversely, growth-independent treatment shifts the mean towards the epithelial phenotype as epithelial cells compensate for the mortality by faster growth. types. For our parameters and assumptions on the growth rates of the phenotypes (Table 1), 285 the switch occurs atλ ≈ −0.089. Consequently, in the one-block adaptive treatment (Fig. 6c), 286 the reductions of more epithelial tumors with λ <λ are identical to reductions obtained for 287 the growth-dependent treatment type. More mesenchymal tumors with λ >λ are aected by 288 one-block adaptive treatment identically to the growth-independent treatment. 289 The eect of sequential multi-block treatment schemes depends on phenotype 290 composition and transition dynamics. 291 To explore the potential of leveraging the treatment-induced phenotype changes, we also 292 investigated sequential treatment patterns that consist of multiple treatment blocks. In the 293 multi-block schemes (see Fig. S2 for tumor burden and mean phenotype dynamics of the two-294 block scheme), the application duration of the optimal treatment is reduced, thus, resulting in a 295 less eective treatment of tumors with a high transition bias (Fig. 6, panels d and e).  For the two-block adaptive treatment, the tumor reduction pattern is more intricate (Fig. 6f ). 301 At λ >λ, the mortality exerted by growth-independent treatment is higher; thus, growth-302 independent treatment is chosen as the rst treatment type. At large transition speed c and 303 high transition bias λ, growth-independent treatment is also chosen for the second treatment 304 block. Only for a transition bias slightly exceedingλ and small c, a treatment switch occurs, 305 and the tumor burden reduction is higher than in the single-block scheme. For λ <λ, the rst 306 treatment type is growth-dependent treatment. Only for small c and close toλ the treatment type 307 is switched after half the treatment period, which interestingly worsens the treatment outcome can aect the predicted reduction in tumor burden. We nd that the tumor burden is reduced 327 more strongly when a single-block growth-dependent treatment is accompanied by an adjuvant 328 drug that makes the tumor more epithelial by decreasing the transition bias λ (arrow in Fig. 6a). 329 Contrastingly, the eect of the single-block growth-independent treatment can be enhanced by 330 an adjuvant drug that makes the tumor more mesenchymal by increasing the transition bias 331 λ (arrow in Fig. 6b). For multi-block treatment schemes, we nd that increasing transition Thus, we nd that phenotypic plasticity can be responsible for both the epithelial-mesenchymal 358 transition at the primary site and the mesenchymal-epithelial transition at the secondary site. 359 However, microenvironmental dierences undoubtedly exacerbate phenotypic changes during the 360 dissemination from primary to secondary site 8 . Indeed, we can capture much more pronounced 361 phenotype changes by assuming dierent transition biases λ 1 ̸ = λ 2 at primary and secondary 362 sites, resulting in diering stable phenotype distributions. 363 Moments of the phenotype distribution guide the choice of phenotype-matched 364 treatment. 365 We investigated growth-dependent and growth-independent treatment types, capturing some 366 mechanisms of chemo-and immunotherapy. We found that the transition bias, and thus the Further, we assume linearly decreasing growth rates from epithelial to mesenchymal 482 phenotypes, reecting the higher proliferation rates of epithelial cells. We scale all growth rates 483 relative to the epithelial growth rate r 1 , and express the growth rate for the i th phenotype as Here, we show the 2 nd (variance) and third central moments of the stable phenotype distribution. For measurements, label the phenotype as i = 1, 2, . . . , N , measure their abundances x i . Mean is, then, the phenotypes i weighted by the measured fractions or relative abundance x i X . a The mean allows the identication of phenotypematched treatment, e.g., to select growth-dependent treatment for mainly epithelial tumors. b The variance represents phenotypic heterogeneity, which often impedes treatment success. c The third central moment is a proxy for the abundance of phenotypes far from the mean phenotype that are potentially unmatched by the phenotype-matched treatment. This measure can thus signal when the mean is an inaccurate representation of the phenotype distribution, for example, when it is strongly skewed. Thus, to maximize the ecacy of the phenotype-matched treatment it is favorable to minimize the heterogeneity and inaccuracy of the mean. High transition bias, and mean phenotype (second row) during and after treatment for tumors with specic transition bias λ (line color) and transition speed c (line type). The growth-dependent treatment block is depicted by a violet bar on the top, and the growth-independent block is shown by a green bar.
The rst block is applied between t = 0 and t = 5, and the second block is applied between t = 5 and t = 10. After the treatment duration, the tumor regrowth is tracked until t = 1000. For the adaptive treatment scheme, the treatment type is chosen at the beginning of the block based on which type exerts the higher mortality on the present phenotype distribution. In practice, to nd the best treatment type, we compare the mortalities exerted by growth-dependent and growth-independent treatment N i=1 m I x i . Their dierence decides the best treatment type, i.e., N i=1 m D r i x i − m I x i > 0 implies that the growth-dependent treatment has higher mortality on the present phenotype distribution. Growth-dependent treatment is chosen for both blocks if the transition bias is low (λ = −0.9, blue lines). For unbiased, slow transitions (λ = 0 and c = 0.01), the treatment type switches from growth-independent to growth-dependent between the rst and the second block (solid gray line). The growth-independent treatment type is applied in both treatment blocks for all other parameter combinations. Compared to single-block treatments (Fig. 5), tumor burden reductions at the end of the treatment phase can be higher or lower in two-block treatment schemes, depending on the transition bias and transition speeds, mirroring the complex patterns in Fig. 6. Two-block treatments are thus not always benecial but can also reduce treatment ecacy if the treatment is not matched with the phenotype distribution.
5.2 Stability of the coexistence equilibrium 806 We investigate the stability of the coexistence phenotype distribution (Eq. (1)) analytically by 807 linear stability analysis for N = 2 and N = 3 57 . For N > 3, we compute the eigenvalues using The following dierential equations describe the system with two phenotypes without treatment.
The matrix allows us to write the system in Eq. S6 asẋ = B(t)x. Since all o-diagonal entries of B(t) 878 are non-negative for t ≥ 0, the system in Eq. S6 is positive (Theorem 3.1) 67 and x(t) ≥ 0 and 879 x(t) ̸ = 0 for all t ≥ 0 meaning that, when the system starts with non-negative initial conditions, 880 it stays non-negative throughout. 881 To study the dynamics of this system, dene the scalar function V as and note that V > 0 for X ̸ = K while V = 0 for X = K. Dierentiating, where from Eq. 2, the denition of X, and telescoping, so that we can write the right-hand side of Eq. S11 in closed form as (S13) This expression shows that the derivative of V is negative whenever X ̸ = K, while it vanishes at 886 X = K. Therefore, V is a Lyapunov function for the dynamics of X and X * = K is a globally 887 stable equilibrium of X. Notably, this stability result is independent of the dynamics of the 888 system that originate from the second term on the right-hand side of Eq. S6. 889 Due to this result, the system in Eq. S6 eventually approaches the linear system 890ẋ = Ax, (S14) as t → ∞, where A is time independent. The system of Eq.
(2) has xed points when the 891 two terms, growth dynamics and transition dynamics, in Eq. (S6) vanish simultaneously. All 892 phenotype conguration satisfying X = K make the rst term vanish. 893 The transition matrix A can be seen as an irreducible generator of a continuous-time Markov  Figure S4: Unequal competition allows tumor burdens higher than the carrying capacity. The panels show approach to the stable phenotype distribution from the same initial condition, (x 1 , x 2 , x 3 ) = ( 1 10 , 0, 0), for dierent combinations of transition speed c and transition bias λ. We assume that cells consume resources proportional to their growth rate, α i = r i · (1 time unit). Since mesenchymal cells consume fewer resources, the higher the mesenchymal proportion, the more cells can be sustained above the carrying capacity. In the original model, we assumed that all phenotypes share the resources identically. However, 919 resource utilization of phenotypes may dier and lead to a dierent competitiveness of 920 phenotypes. To account for this aspect, here we extend the model by including phenotype-921 specic competition coecients α i . Thus, we replace X in Eqs. (2) withX = N j=1 α j x j and 922 obtain for the untreated dynamics The intensity of growth-independent treatment is m I ≈ 0.65. Darker colors represent a higher reduction and, thus, a better outcome.
We evaluate the eect of splitting the treatment period into multiple treatment blocks (rows) and investigate dierent treatment schemes with either predened or adaptive treatment sequences (columns). distribution changes with the transition bias λ. When there is no transition bias to switch to either an epithelial or a mesenchymal-like phenotype, i.e., λ = 0, the stable distribution is not uniform as faster proliferating cells transition faster. Transition bias towards the epithelial-like phenotype, λ < 0, leads to an increase in epithelial cells. Conversely, a transition bias towards mesenchymal-like phenotypes, λ > 0, leads to an increase in mesenchymal cells. Notably, now the abundances, particularly of the mesenchymal phenotype, can exceed the 935 carrying capacity K, i.e. κ > 1, up to an eective carrying capacity given byK = κK 936 (Fig. S4). Consider the case where mesenchymal-like cells use fewer resources than epithelial-like 937 cells. Namely, we assume that the competition coecients are proportional to the growth rate. 938 Rescaling the competition coecients then gives α i = r i · 1 (time unit) for the i th phenotype. 939 Note that the units of growth rate and competition coecients are dierent requiring the 940 multiplication by 1 (time unit). This assumption incorporates a classical r − K trade-o into 941 our model with more proliferative but less competitive epithelial cells and less proliferative but 942 more competitive mesenchymal cells 69 . This higher abundance of less proliferative but more 943 competitive mesenchymal phenotypes increases the treatment eciencies for intermediate and 944 high transition biases (Fig. S5). Similarly, when a phenotype is harmful to others α j > 1 the 945 equilibrium abundance can be lower than the carrying capacity κ < 1.
(S16) 974 This set of equations has two equilibria similar to the original model. One equilibrium 975 is a saddle-node equilibrium of extinction, and the other is a stable equilibrium with a xed 976 phenotype distribution. This stable equilibrium is given byx * (1−λ) N −j (1+λ) j−1 .

977
The model dynamics with phenotype-dependent transition speed are qualitatively similar to 978 the original model. In this case, however, the stable phenotype distribution depends on the 979 phenotype growth rate, increasing the fraction of more mesenchymal phenotypes (Fig. S6). This 980 leads to stable phenotype distributions with negative third central moment for a larger transition 981 bias range (Fig. S7) and an eciency loss of purely growth-dependent treatment (Fig. S8).